% m_eval_1945_50.m


clear all; 
close all;
format shortG;

% ==================== LOAD FUNDAMENTAL AND DATA  =========================
load parameters.mat;                     % basic parameters 
load calib_productivity_amenity.mat;     % inverted fundamentals every 5 years from 1955
load calib_setup.mat;                    % Info data every 5 years from 1945
load commuting_cost.mat;  
load calib_auxiliary.mat;                % Info used in this calibration 

Rdata45 = Rdata_all(:,1); Rdata50 = Rdata_all(:,2); 
Ldata45 = Ldata_all(:,1); Ldata50 = Ldata_all(:,2); 
N = size(Ldata_all,1); 

% ==================== COMPUTE COMMUTING PATTERN 1945 =====================
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-5,'TolX',1e-5);
K45=Kappa_mat(:,:,1).^(-rho/sigma); 
Y0=ones(N,1);     
comsyst = @(y)fcomeval2(y,K45,Rdata45,Ldata45,rho,sigma);
[y_fsolve, yfval_fsolve]=fsolve(comsyst, Y0, options0); Y45=abs(y_fsolve); Y45=Y45./geomean(Y45);  
ytemp=K45.*(repmat(Y45,1,N).^(rho/sigma)); 
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
Lcom45=ytemp.*repmat(Rdata45',N,1); 

% =================== PARAMETERS ==========================================
alpha=alphaest;   beta=betaest; 

% =========== COMPUTE LOWER BOUNDS for Theta 1950 & SET Theta 1950 ======== 
Rthetalb = (Rdata50 - Rdata45 < 0);   Rthetalb = (1 - Rdata50./Rdata45).*Rthetalb;  Rthetalb = max(Rthetalb); 
Lthetalb = (Ldata50 - Ldata45 < 0);   Lthetalb = (1 - Ldata50./Ldata45).*Lthetalb;  Lthetalb = max(Lthetalb);
thetalb = ceil(100.*max([Rthetalb; Lthetalb]))./100;

theta50 = 0.90;
if min(theta50 - thetalb)<0
    disp('>>> CHECK THETA VALUES (LOWER BOUND) ! <<< ');   pause
end 
theta55 = thetavec(1); 

%  ================ INVERSION OF (X,Y) FOR 1950 =============== 
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-6,'TolX',1e-6);
X0=ones(N,1);  
Y0=ones(N,1);  
Rsyst = @(x)fReval2(x,K50,Q50,Rdata45,Rdata50,Ldata45,Ldata50,rho50,sigma50,theta50,mu,gammaL,gammaH);
Lsyst = @(y)fLeval2(y,K50,Q50,Rdata45,Rdata50,Ldata45,Ldata50,rho50,sigma50,theta50,mu,gammaL,gammaH);  
[x_fsolve, xfval_fsolve]=fsolve(Rsyst, X0, options0);   
[y_fsolve, yfval_fsolve]=fsolve(Lsyst, Y0, options0); 
X50=abs(x_fsolve); X50=X50./geomean(X50);
Y50=abs(y_fsolve); Y50=Y50./geomean(Y50);

% ************************************************
% *********** FORWARD LOOKING MODEL **************
% ************************************************
% =========== DECOMPOSE (X,Y) INTO FUNDAMENTALS AND AGGLOMERATION =========
% compute fundamental productivity and amenities in 1950
b50 = X50.*((Rdata50./Area).^(-beta)).*(X55.^(-rho*(1-theta55))); 
a50 = Y50.*((Ldata50./Area).^(-alpha)).*(Y55.^(-rho*(1-theta55))); 

% ========== COMPUTE STRUCTURAL RESIDUALS (FORWARD LOOKING) ================
l_ave_b = sum(log(bb),2)./size(bb,2); % average amenities over time (1955-75) 
l_ave_a = sum(log(aa),2)./size(aa,2); % average productivity over time (1955-75)

l_bb50 = log(b50) - l_ave_b;         % (log) structural error in amenities in 1950 (Forward Looking)
l_aa50 = log(a50) - l_ave_a;         % (log) structural error in productivity in 1950 (Forward Looking)

% Note: Since we have included the time-fixed fundamental for 1950, this
% is actually fundamentals + 1950 year fixed effects. 
l_b50_ave = sum(log(b50))./N; 
l_a50_ave = sum(log(a50))./N; 
bb_err50 = exp(log(b50)-l_b50_ave); 
aa_err50 = exp(log(a50)-l_a50_ave); 

XX50 = exp(l_ave_b).*((Rdata50./Area).^beta).*(X55.^(rho*(1-theta55)));   
YY50 = exp(l_ave_a).*((Ldata50./Area).^alpha).*(Y55.^(rho*(1-theta55))); 

% ************************************************
% ******* FIGURE FOR STRUCTURAL RESIDUALS ********
% ***** CHECK MOMENT CONDITIONS (FIGURE 5A) ******
% ************************************************
% Store residuals
d_aa_err = log(aa_err(:,2:5))-log(aa_err(:,1:4)); 
d_bb_err = log(bb_err(:,2:5))-log(bb_err(:,1:4)); 

d_aa_err_mean = sum(d_aa_err,2)./size(d_aa_err,2); 
d_bb_err_mean = sum(d_bb_err,2)./size(d_bb_err,2);

mm_aa_res = Gmat*d_aa_err_mean./sum(Gmat,2); 
mm_bb_res = Gmat*d_bb_err_mean./sum(Gmat,2); 
mm_aa_res = Gmat'*mm_aa_res; 
mm_bb_res = Gmat'*mm_bb_res; 

res_difmat = [d_aa_err_mean mm_aa_res d_bb_err_mean mm_bb_res cbddist]; 
res_difmat2 = sortrows(res_difmat, 5); 

% Figure based on Rank (x-axis)
figure;
scatter((1:N)', res_difmat2(:,1), 14, 'ok','filled');
ylabel('Change in log of residual fundamentals', 'fontsize',12);
yticks([-0.2 -0.1 0 0.1 0.2]);
set(gca,'YColor','k');
hold on;  %yyaxis right
plot((1:N)', res_difmat2(:,2), '-', 'Color', '#9e003a', 'LineWidth', 3.0); 
hold on;
scatter((1:N)', res_difmat2(:,3), 14, [0.7 0.7 0.7], 'filled', 'square');
set(gca,'YColor','k');
hold on;  %yyaxis right
plot((1:N)', res_difmat2(:,4), '-', 'Color', '#4DBEEE', 'LineWidth', 2.0); 
xlim([0 N+1]); ylim([-0.2 0.2]); hold off;
xlabel('Locations Ranked by CBD distance (from CBD to Suburb)','fontsize',12); %from short to long 
legend('Fundamental productivity changes','Average of fundamental productivity changes across grids','Fundamental amenity changes','Average of fundamental amenity changes across grids');
set(gca,'YColor','k');
f=gca;
exportgraphics(f,'../../output/figure/figure_momentcondition_check.pdf');
close all; 

% ***********************************************************
% ******* FIGURE FOR STRUCTURAL RESIDUALS IN 1950-55 ********
% ***********************************************************
d_aa_err_50 = log(aa_err(:,1))-log(aa_err50);
d_bb_err_50 = log(bb_err(:,1))-log(bb_err50);
res_difmat_50 = [d_aa_err_50 d_bb_err_50 cbddist]; 
res_difmat_50 = sortrows(res_difmat_50, 3); 

figure;
scatter((1:N)', res_difmat_50(:,1), 15, 'ok','filled');
ylabel('Change in log of residual fundamentals in 1950-55', 'fontsize',12);
yticks([-0.2 -0.1 0 0.1 0.2]);
set(gca,'YColor','k'); hold on 
scatter((1:N)', res_difmat_50(:,2), 20, [0.7 0.7 0.7], 'filled', 'square');
set(gca,'YColor','k');
xlim([0 N+1]); ylim([-0.2 0.2]); hold off;
xlabel('Locations Ranked by CBD distance (from CBD to Suburb)','fontsize',12); %from Low density to High density
legend('Productivity 1950-55','Amenities 1950-55'); 
set(gca,'YColor','k');
f=gca;
exportgraphics(f,'../../output/figure/figure_fundamental_momentcondition_50_55.pdf');
close all; 

% ********************************************************************
% ******** SIMULATE FORWARD LOOKING MODEL ****************************
% ******** WITHOUT STRUCTURAL RESIDUALS ******************************
% ** NOTE: HERE WE USE PROBABILITIES THAT IS BASED ON XX50, YY50 *****
% ********************************************************************
XX_temp = K50.*((repmat(Q50',N,1).^(-mu)).*repmat(XX50',N,1)).^(rho50/sigma50); 
XX_temp = XX_temp./repmat(sum(XX_temp,2),1,N); 

YY_temp = K50.*((repmat(Q50,1,N).^(-gammaH/gammaL)).*(repmat(YY50,1,N).^(1/gammaL))).^(rho50/sigma50); 
YY_temp = YY_temp./repmat(sum(YY_temp,1),N,1);

R50_p=Rdata50; 
L50_p=Ldata50; 
maxit=1e+6; 
ii=0; 
update=0.05;
while ii < maxit
    Ldif_p = L50_p-(1-theta50).*Ldata45;
    Rdif_p = R50_p-(1-theta50).*Rdata45; 
    
    R50_n = (1-theta50).*Rdata45+sum(XX_temp.*repmat(Ldif_p,1,N),1)'; 
    L50_n = (1-theta50).*Ldata45+sum(YY_temp.*repmat(Rdif_p',N,1),2);
    R50_p_r = round(R50_p*1e+6)./1e+6; R50_n_r = round(R50_n*1e+6)./1e+6; 
    L50_p_r = round(L50_p*1e+6)./1e+6; L50_n_r = round(L50_n*1e+6)./1e+6; 

    if min(R50_p_r == R50_n_r)==1 && min(L50_p_r == L50_n_r)==1
        break; 
    else
        ii=ii+1; 
        R50_p = (1-update).*R50_p+update.*R50_n; 
        L50_p = (1-update).*L50_p+update.*L50_n;
        if rem(ii,100)==1
        disp(['Largest Error is:',num2str(ii), ' ',num2str(max(abs(R50_p_r-R50_n_r))), ' ',num2str(max(abs(L50_p_r-L50_n_r)))]);
        end
    end
end 

% ============ WELFARE COMPARISON WITH DATA ===============================
% compute commuting patterns in 1950
LQ50 = LQmat(:,1); 
lambda50 = K50.*((repmat(LQ50',N,1).^(-mu)).*repmat(X50',N,1)).^(rho50/sigma50); 
lambda50 = lambda50./repmat(sum(lambda50,2),1,N); 
Lcom50 = (1-theta50).*Lcom45+lambda50.*repmat(Ldata50-(1-theta50).*Ldata45,1,N); 
V50=K50.*((repmat(LQ50',N,1).^(-mu)).*repmat(X50',N,1)).*((repmat(LQ50,1,N).^(-gammaH/gammaL)).*(repmat(Y50,1,N).^(1/gammaL)));  % V_{in,50}
VV50=sigma50*log(sum(sum((Lcom50./sum(sum(Lcom50))).*(exp(V50).^(rho50/sigma50))))); % V_t eq.(G.37)

lambda50_p=K50.*((repmat(LQ50',N,1).^(-mu)).*repmat(XX50',N,1)).^(rho50/sigma50); 
lambda50_p = lambda50_p./repmat(sum(lambda50_p,2),1,N); 
Lcom50_p = (1-theta50).*Lcom45+lambda50_p.*repmat(L50_p-(1-theta50).*Ldata45,1,N); 
V50_p=K50.*((repmat(LQ50',N,1).^(-mu)).*repmat(XX50',N,1)).*((repmat(LQ50,1,N).^(-gammaH/gammaL)).*(repmat(YY50,1,N).^(1/gammaL)));
VV50_p=sigma50*log(sum(sum((Lcom50_p./sum(sum(Lcom50_p))).*(exp(V50_p).^(rho50/sigma50)))));

VV50_dif=VV50_p-VV50; 
VV50_dif_percent=VV50_dif/rho50; 

% ============ VARIANCE DECOMPOSITION FOR FORWARD LOOKING =================
v_X50 = var(log(X50)); 
v_Y50 = var(log(Y50));
cov_X50 = cov(log(XX50),l_bb50); 
cov_Y50 = cov(log(YY50),l_aa50); 

% ==================  SAVE RESULTS ========================================
disp(' >>> ====== FINISH 1945-50 CALIBRATION AND SAVE RESULTS ====== <<< ')

varresult=[v_X50, v_Y50, cov_X50(1,1),cov_Y50(1,1),cov_X50(2,2),cov_Y50(2,2),cov_X50(1,2), cov_Y50(1,2)];
resultparam=array2table(varresult,'VariableNames',{'var_X','var_Y','var_XX50', ... 
    'var_YY50','var_aa50','var_bb50', 'cov_X50','cov_Y50'}); 
writetable(resultparam,'../../output/quant/table_fundamentals_variance_1950.csv'); 

result=array2table([ID, X50, Y50, XX50, YY50, Rdata45, Rdata50, Ldata45,Ldata50, ...
    R50_p, L50_p, cbddist, Area],'VariableNames',{'geocode','X50','Y50','XX50','YY50',...
    'pop45','pop50','emp45','emp50','pop50_pred','emp50_pred','cbddist', 'area'}); 
writetable(result,'../../output/quant/output_evaluations_1945_50.csv');

save('fundamentals_1950.mat','b50','a50','theta50'); 

fundamentalresult=array2table([a50 b50 aa_err50 bb_err50  cbddist, Area, ID, lat, lon],'VariableNames', {'aa50','bb50',... 
    'aa_err50','bb_err50','cbddist', 'area', 'geocode', 'lat', 'lon'}); 
writetable(fundamentalresult,'../../output/quant/output_fundamentals_1945_50.csv');

welfareresult=array2table([VV50 VV50_p VV50_dif VV50_dif_percent],'VariableNames', {'avewelfare50','avewelfare50_model',... 
    'avewelfare50_dif', 'avewelfare50_percentagedif'}); 
writetable(welfareresult,'../../output/quant/output_avewelfare_1945_50.csv');

save('calib_1950_auxiliary.mat','X50','Y50','a50','b50'); 

% ************************************************************************
% ************************************************************************
% **************************  FUNCTIONS **********************************
% ************************************************************************
% ************************************************************************
function[Reval] = fReval2(X,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH)   % eq D.16. X corresponds to Theta

N = size(DRpre,1); 
X=abs(X);  % avoid negative values
X=X./geomean(X);
xtemp=K.*((repmat(Q',N,1).^(-mu)).*repmat(X',N,1)).^(rho/sigma);
xtemp=xtemp./repmat(sum(xtemp,2),1,N); 
Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 
Reval = Rdif - sum(xtemp.*repmat(Ldif,1,N),1)'; 
Reval = abs(Reval); 

end

function[Leval] = fLeval2(Y,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH) % eq D.17. Y corersponds to Omega

N = size(DRpre,1); 
Y=abs(Y);  % avoid negative values
Y=Y./geomean(Y);
ytemp= K.*((repmat(Q,1,N).^(-gammaH/gammaL)).*(repmat(Y,1,N).^(1/gammaL))).^(rho/sigma);
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 
Leval = Ldif - sum(ytemp.*repmat(Rdif',N,1),2); 
Leval = abs(Leval); 

end

function[comeval] = fcomeval2(Y,K,DR,DL,rho,sigma) % eq D.17. Y corersponds to Omega

N = size(DR,1); 
Y=abs(Y);  % avoid negative values
Y=Y./geomean(Y);
ytemp=K.*(repmat(Y,1,N).^(rho/sigma)); 
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
comeval = DL - sum(ytemp.*repmat(DR',N,1),2); 
comeval = abs(comeval); 

end
